package fem2;

import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.IMatrix;

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

import math.GaussPoints;
import math2.GaussTetrahedra;
import math2.GaussTriangleGid;
import math2.LinAlgUtilities;
import math2.delaunay.Triangle;
import fem2.element.StructuralElement;
import fem2.enu.IntegrationType;
import fem2.enu.MeshPartType;

/**
 * singleton class to manage all intergration points for all type of geometry
 * 
 * @author hbui
 * 
 */
public class IntegrationManager {

	private static IntegrationManager instance = new IntegrationManager();

	private IntegrationManager() {
	}

	/**
	 * Get the instance of IntegrationManager singleton
	 * 
	 * @return the Singleton instance
	 */
	public static IntegrationManager getInstance() {
		return instance;
	}

	/**
	 * Utility method to set all elements in the model to a specific integration
	 * type
	 * 
	 * @param M
	 * @param type
	 */
	public void setIntegrationType(Model M, IntegrationType type) {
		for (int i = 0; i < M.countElements(); i++) {
			Element e = M.getElement(i);
			if (e instanceof StructuralElement) {
				((StructuralElement) e).setIntegrationType(type);
			}
		}
	}

	/**
	 * check if the xi fall in the reference domain of mesh part which is
	 * postulated by integration rule
	 * 
	 * @param type
	 * @param xi
	 * @return
	 */
	public boolean checkBoundary(MeshPartType type, double... xi) {
		if (type == MeshPartType.TRI3 || type == MeshPartType.TRI6) {
			/*
			 * for triangle the boundary of reference domain is {(x,y) |
			 * 0<=x,y,x+y<=1
			 */
			return (0 <= xi[0]) && (xi[0] <= 1) && (0 <= xi[1]) && (xi[1] <= 1)
					&& (0 <= xi[0] + xi[1]) && (xi[0] + xi[1] <= 1);
		} else if (type == MeshPartType.QUAD4 || type == MeshPartType.QUAD8
				|| type == MeshPartType.QUAD9) {
			/*
			 * for quadrilateral the boundary of reference domain is {(x,y) |
			 * -1<=x,y<=1
			 */
			return (-1 <= xi[0]) && (xi[0] <= 1) && (-1 <= xi[1]) && (xi[1] <= 1);
		} else if (type == MeshPartType.TET4 || type == MeshPartType.TET10) {
			/*
			 * for tetrahedra the boundary of reference domain is {(x,y,z) |
			 * 0<=x,y,z,x+y+z<=1
			 */
			return (0 <= xi[0]) && (xi[0] <= 1) && (0 <= xi[1]) && (xi[1] <= 1) && (0 <= xi[2])
					&& (xi[2] <= 1) && (0 <= xi[0] + xi[1] + xi[2]) && (xi[0] + xi[1] + xi[2] <= 1);
		} else if (type == MeshPartType.HEX8 || type == MeshPartType.HEX20) {
			/*
			 * for hexahedra the boundary of reference domain is {(x,y,z) |
			 * -1<=x,y,z<=1
			 */
			return (-1 <= xi[0]) && (xi[0] <= 1) && (-1 <= xi[1]) && (xi[1] <= 1) && (-1 <= xi[2])
					&& (xi[2] <= 1);
		} else {
			throw new Error("this mesh part type is not supported");
		}
	}

	public double[] getReferenceCoordinate(MeshPartType type, int i) {
		switch (type) {
		case TRI3:
			switch (i) {
			case 0:
				return new double[] { 0.0, 0.0 };
			case 1:
				return new double[] { 1.0, 0.0 };
			case 2:
				return new double[] { 0.0, 1.0 };
			default:
				throw new Error("index exceeds number of nodes on the mesh part");
			}
		case QUAD4:
			switch (i) {
			case 0:
				return new double[] { -1.0, -1.0 };
			case 1:
				return new double[] { 1.0, -1.0 };
			case 2:
				return new double[] { 1.0, 1.0 };
			case 3:
				return new double[] { -1.0, 1.0 };
			default:
				throw new Error("index exceeds number of nodes on the mesh part");
			}
			/*
			 * TODO: extends to other mesh part type
			 */
		default:
			throw new Error("this mesh part type is not supported");
		}
	}

	/***********************************************
	 * Utility method for selective reduced integration scheme
	 ***********************************************/

	double[][] D2 = new double[][] { //
	{ 1.0 / 2, 1.0 / 2, 0.0 }, //
			{ 1.0 / 2, 1.0 / 2, 0.0 }, //
			{ 0.0, 0.0, 0.0 } };
	double[][] ID2 = new double[][] { //
	{ 1.0 / 2, -1.0 / 2, 0.0 }, //
			{ -1.0 / 2, 1.0 / 2, 0.0 }, //
			{ 0.0, 0.0, 1.0 } };
	double[][] D3 = new double[][] { { 1.0 / 3, 1.0 / 3, 1.0 / 3, 0.0, 0.0, 0.0 },
			{ 1.0 / 3, 1.0 / 3, 1.0 / 3, 0.0, 0.0, 0.0 },
			{ 1.0 / 3, 1.0 / 3, 1.0 / 3, 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, //
			{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, //
			{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } };
	double[][] ID3 = new double[][] { { 2.0 / 3, -1.0 / 3, -1.0 / 3, 0.0, 0.0, 0.0 },
			{ -1.0 / 3, 2.0 / 3, -1.0 / 3, 0.0, 0.0, 0.0 },
			{ -1.0 / 3, -1.0 / 3, 2.0 / 3, 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 }, //
			{ 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 }, //
			{ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };

	private IMatrix volumIntegralDMatrix2D = new Array2DMatrix(D2);
	private IMatrix volumIntegralIminusDMatrix2D = new Array2DMatrix(ID2);
	private IMatrix volumIntegralDMatrix3D = new Array2DMatrix(D3);
	private IMatrix volumIntegralIminusDMatrix3D = new Array2DMatrix(ID3);

	public IMatrix getVolumeIntegralDMatrix(int dim) {
		if (dim == 2)
			return volumIntegralDMatrix2D;
		else if (dim == 3)
			return volumIntegralDMatrix3D;
		else {
			throw new Error("invalid dimension");
		}
	}

	public IMatrix getVolumeIntegralIminusDMatrix(int dim) {
		if (dim == 2)
			return volumIntegralIminusDMatrix2D;
		else if (dim == 3)
			return volumIntegralIminusDMatrix3D;
		else {
			throw new Error("invalid dimension");
		}
	}

	/***********************************************
	 * trivial integration points provider on mesh part
	 ***********************************************/

	public double[] getIntegrationPoint(MeshPartType mpType, int p, int k) {
		if (mpType == MeshPartType.LINE2 || mpType == MeshPartType.LINE3) {
			return new double[] { GaussPoints.getPoint(p, k) };
		} else if (mpType == MeshPartType.TRI3 || mpType == MeshPartType.TRI6) {
			// return new double[] {
			// (1 + GaussTrianglePoints.getGaussXi1(p, k)) / 2,
			// (1 + GaussTrianglePoints.getGaussXi2(p, k)) / 2 };
			return new double[] { GaussTriangleGid.getGaussXi1(p, k),
					GaussTriangleGid.getGaussXi2(p, k) };
		} else if (mpType == MeshPartType.QUAD4 || mpType == MeshPartType.QUAD8) {
			int i = k / p;
			int j = k % p;
			return new double[] { GaussPoints.getPoint(p, i), GaussPoints.getPoint(p, j) };
		} else if (mpType == MeshPartType.TET4 || mpType == MeshPartType.TET10) {
			return new double[] { GaussTetrahedra.getGaussXi1(p, k),
					GaussTetrahedra.getGaussXi2(p, k), GaussTetrahedra.getGaussXi3(p, k) };
		} else if (mpType == MeshPartType.HEX8 || mpType == MeshPartType.HEX20) {
			int l = k % p;
			int j = (k / p) % p;
			int i = (k / p) / p;
			return new double[] { GaussPoints.getPoint(p, i), GaussPoints.getPoint(p, j),
					GaussPoints.getPoint(p, l) };
		} else {
			throw new Error("this mesh part type is not supported");
		}
	}

	public double getIntegrationWeight(MeshPartType mpType, int p, int k) {
		if (mpType == MeshPartType.LINE2 || mpType == MeshPartType.LINE3) {
			return GaussPoints.getWeight(p, k);
		} else if (mpType == MeshPartType.TRI3 || mpType == MeshPartType.TRI6) {
			// return GaussTrianglePoints.getGaussWeight(p, k) / 4;
			return GaussTriangleGid.getGaussWeight(p, k);
		} else if (mpType == MeshPartType.QUAD4 || mpType == MeshPartType.QUAD8) {
			int i = k / p;
			int j = k % p;
			double w1 = GaussPoints.getWeight(p, i);
			double w2 = GaussPoints.getWeight(p, j);
			return w1 * w2;
		} else if (mpType == MeshPartType.TET4 || mpType == MeshPartType.TET10) {
			return GaussTetrahedra.getGaussWeight(p, k);
		} else if (mpType == MeshPartType.HEX8 || mpType == MeshPartType.HEX20) {
			int l = k % p;
			int j = (k / p) % p;
			int i = (k / p) / p;
			double w1 = GaussPoints.getWeight(p, i);
			double w2 = GaussPoints.getWeight(p, j);
			double w3 = GaussPoints.getWeight(p, l);
			return w1 * w2 * w3;
		} else {
			throw new Error("this mesh part type is not supported");
		}
	}

	public int getNumIntegrationPoints(MeshPartType mpType, int p) {
		if (mpType == MeshPartType.LINE2 || mpType == MeshPartType.LINE3) {
			return p;
		} else if (mpType == MeshPartType.TRI3 || mpType == MeshPartType.TRI6) {
			// return GaussTrianglePoints.getNumGaussPoints(p);
			return GaussTriangleGid.getNumGaussPoints(p);
		} else if (mpType == MeshPartType.QUAD4 || mpType == MeshPartType.QUAD8) {
			return p * p;
		} else if (mpType == MeshPartType.TET4 || mpType == MeshPartType.TET10) {
			return GaussTetrahedra.getNumGaussPoints(p);
		} else if (mpType == MeshPartType.HEX8 || mpType == MeshPartType.HEX20) {
			return p * p * p;
		} else {
			throw new Error("this mesh part type is not supported");
		}
	}

	/***********************************************
	 * integration points for xfem
	 ***********************************************/

	private Map<Element, ArrayList<IntegrationPoint>> elementIntegrationMap = new HashMap<Element, ArrayList<IntegrationPoint>>();

	/**
	 * associate element with list of integration points
	 * 
	 * @param e
	 * @param pointList
	 */
	public void setIntegrationPoints(Element e, ArrayList<IntegrationPoint> pointList) {
		elementIntegrationMap.put(e, pointList);
	}

	/**
	 * get number of integration points on element
	 * 
	 * @param e
	 * @return
	 */
	public int getNumIntegrationPoints(Element e) {
		return elementIntegrationMap.get(e).size();
	}

	/**
	 * get integration point
	 * 
	 * @param e
	 * @return
	 */
	public double[] getIntegrationPoint(Element e, int k) {
		return elementIntegrationMap.get(e).get(k).getPoint();
	}

	/**
	 * get integration weight
	 * 
	 * @param e
	 * @return
	 */
	public double getIntegrationWeight(Element e, int k) {
		return elementIntegrationMap.get(e).get(k).getWeight();
	}

	/**
	 * create integration points for sub-triangle mesh
	 * 
	 * @param order
	 * @param triangles
	 * @return
	 */
	public ArrayList<IntegrationPoint> createIntegrationPoints(int order, List<Triangle> triangles) {
		IntegrationManager im = IntegrationManager.getInstance();
		ArrayList<IntegrationPoint> ips = new ArrayList<IntegrationPoint>();
		for (Triangle t : triangles) {
			Face f;
			double[] xi1 = t.getPi().getCoordinate();
			double[] xi2 = t.getPj().getCoordinate();
			double[] xi3 = t.getPk().getCoordinate();
			// if (order == 1) {
			Node n1 = new Node(xi1);
			Node n2 = new Node(xi2);
			Node n3 = new Node(xi3);
			f = new Face(n1, n2, n3);
			// } else if (order == 2) {
			// Node n1 = new Node(xi1);
			// Node n2 = new Node(xi2);
			// Node n3 = new Node(xi3);
			// Node n4 = new Node(0.5 * (xi1[0] + xi2[0]), 0.5 * (xi1[1] +
			// xi2[1]));
			// Node n5 = new Node(0.5 * (xi2[0] + xi3[0]), 0.5 * (xi2[1] +
			// xi3[1]));
			// Node n6 = new Node(0.5 * (xi3[0] + xi1[0]), 0.5 * (xi3[1] +
			// xi1[1]));
			// f = new Face(n1, n2, n3, n4, n5, n6);
			// } else {
			// throw new Error("higher order integration is not supported");
			// }
			int ng = im.getNumIntegrationPoints(f.getType(), order);
			for (int k = 0; k < ng; k++) {
				double[] xiLocal = im.getIntegrationPoint(f.getType(), order, k);
				double[] ip = f.XAt(xiLocal);
				double detJ = LinAlgUtilities.determinantNxN(f.jacXAt(xiLocal));
				double iw = im.getIntegrationWeight(f.getType(), order, k) * detJ;
				ips.add(new IntegrationPoint(ip, iw));
			}
		}
		return ips;
	}

	/**
	 * get default integration order for a mesh part type
	 * 
	 * @param type
	 * @return
	 */
	public int getDefaultIntegrationOrder(MeshPartType type) {
		switch (type) {
		case TRI3:
			return 1;
		case TRI6:
			return 2;
		case QUAD4:
			return 2;
		case QUAD8:
			return 3;
		case QUAD9:
			return 3;
		case TET4:
			return 1;
		case TET10:
			return 2;
		case HEX8:
			return 2;
		case HEX20:
			return 3;
		case LINE2:
			return 1;
		case LINE3:
			return 2;
		default:
			throw new Error("not supported mesh part type");
		}
	}
}
